Brief description of the script
Content
This script reads and plots the data from the 3D and 2D morphometric
analysis of the artefacts’ engraved incisions. For details on the
methods and data acquisition, please visit the Materials and Methods
section of the paper.
The knit directory for this script is the project directory.This R
project and respective scripts follow the procedures described by
Marwick et al. 2017.
For any questions, comments and inputs, please contact:
João Marreiros, joao.marreiros@leiza.de
Load packages
library(R.utils)
## Warning: package 'R.utils' was built under R version 4.3.1
## Loading required package: R.oo
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.2 (2022-06-13 22:00:14 UTC) successfully loaded. See ?R.methodsS3 for help.
## R.oo v1.25.0 (2022-06-12 02:20:02 UTC) successfully loaded. See ?R.oo for help.
##
## Attaching package: 'R.oo'
## The following object is masked from 'package:R.methodsS3':
##
## throw
## The following objects are masked from 'package:methods':
##
## getClasses, getMethods
## The following objects are masked from 'package:base':
##
## attach, detach, load, save
## R.utils v2.12.3 (2023-11-18 01:00:02 UTC) successfully loaded. See ?R.utils for help.
##
## Attaching package: 'R.utils'
## The following object is masked from 'package:utils':
##
## timestamp
## The following objects are masked from 'package:base':
##
## cat, commandArgs, getOption, isOpen, nullfile, parse, warnings
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.1
library(tools)
library(tidyverse)
## Warning: package 'readr' was built under R version 4.3.1
## Warning: package 'dplyr' was built under R version 4.3.1
## Warning: package 'stringr' was built under R version 4.3.1
## Warning: package 'lubridate' was built under R version 4.3.1
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks R.utils::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(doBy)
## Warning: package 'doBy' was built under R version 4.3.1
##
## Attaching package: 'doBy'
##
## The following object is masked from 'package:dplyr':
##
## order_by
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 4.3.1
library(flextable)
## Warning: package 'flextable' was built under R version 4.3.1
##
## Attaching package: 'flextable'
##
## The following object is masked from 'package:purrr':
##
## compose
library(readr)
library(ggpubr)
##
## Attaching package: 'ggpubr'
##
## The following objects are masked from 'package:flextable':
##
## border, font, rotate
library(crosstable)
## Warning: package 'crosstable' was built under R version 4.3.1
##
## Attaching package: 'crosstable'
##
## The following object is masked from 'package:purrr':
##
## compact
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
Import and preview data
db <- read.csv("../raw_data/data.csv", stringsAsFactors=FALSE, sep=";", dec=".")
str(db)
## 'data.frame': 324 obs. of 8 variables:
## $ site : chr "Amud1" "Qafzeh" "Amud1" "Amud1" ...
## $ aoi : int 1 2 1 2 1 1 1 2 1 2 ...
## $ incision : int 11 3 20 8 11 4 20 10 8 9 ...
## $ profile : int 2 2 1 1 1 3 3 2 3 3 ...
## $ horizontaldistance: num 253 973 300 293 261 ...
## $ maximumdepth : num 10.62 9.05 11.61 12.49 14.43 ...
## $ holearea : int 1366 1378 1710 1718 1877 2043 2048 2110 2124 2237 ...
## $ angle : num 163 148 167 162 162 ...
head(db)
## site aoi incision profile horizontaldistance maximumdepth holearea angle
## 1 Amud1 1 11 2 253.2 10.62 1366 163.1
## 2 Qafzeh 2 3 2 973.1 9.05 1378 148.2
## 3 Amud1 1 20 1 300.3 11.61 1710 167.1
## 4 Amud1 2 8 1 293.2 12.49 1718 161.8
## 5 Amud1 1 11 1 261.4 14.43 1877 162.5
## 6 Amud1 1 4 3 277.0 16.01 2043 158.3
Import and summarize data
# summarise data frome ach artefact, organised by AOI
# sample data per artefact
manotdb <- filter(db, site == "Manot")
amud1db <- filter(db, site == "Amud1")
amud2db <- filter(db, site == "Amud2")
qzdb <- filter(db, site =="Qafzeh")
qnrtdb <- filter(db, site =="Quneitra")
# Manot
# Horizontal distance
manotdistance <- manotdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
manotdepth <- manotdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
manotarea <- manotdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
manotangle <- manotdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
manotdistance
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 9 929. 690. 826. 82.1 814.
## 2 2 12 891. 803. 848. 33.3 848.
## 3 3 42 1163. 587. 811. 125. 816.
manotdepth
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 9 76.5 51.0 62.4 9.11 57.4
## 2 2 12 59.1 44.4 52.3 4.14 53.4
## 3 3 42 120. 25.5 64.8 21.1 61.6
manotarea
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 1 9 34240 10820 25592. 7255. 24590
## 2 2 12 33056 9471 20100. 10597. 18449
## 3 3 42 60785 12293 29378. 10877. 26268.
manotangle
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 9 158. 143. 150. 5.27 151.
## 2 2 12 153. 147 151. 1.65 151
## 3 3 42 160. 136 148. 5.95 148.
# save the results
write_csv(manotdepth, "../summary_stats/manotdepth.csv")
write_csv(manotarea, "../summary_stats/manotarea.csv")
write_csv(manotangle, "../summary_stats/manotangle.csv")
write_csv(manotdistance, "../summary_stats/manotdistance.csv")
# Amud1
# Horizontal distance
amud1distance <- amud1db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
amud1depth <- amud1db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
amud1area <- amud1db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
amud1angle <- amud1db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
amud1distance
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 63 718 208. 393. 119. 370
## 2 2 30 432. 241. 343. 56.3 345
## 3 3 18 506. 323. 430. 58.8 440.
amud1depth
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 63 82.5 5.87 35.2 18.0 32.8
## 2 2 30 56.5 12.5 29.5 11.2 27.5
## 3 3 18 208. 25.9 62.1 41.3 51.2
amud1area
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 1 63 30636 549 8115. 6743. 5774
## 2 2 30 11676 1718 5465. 2525. 5014
## 3 3 18 25652 4953 12850. 6061. 12006.
amud1angle
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 63 168. 123. 152. 9.24 154.
## 2 2 30 166. 131. 154. 8.64 156.
## 3 3 18 170. 149. 158. 6.24 159.
# save the results
write_csv(amud1depth, "../summary_stats/amud1depth.csv")
write_csv(amud1area, "../summary_stats/amud1area.csv")
write_csv(amud1angle, "../summary_stats/amud1angle.csv")
write_csv(amud1distance, "../summary_stats/amud1distance.csv")
# Amud2
# Horizontal distance
amud2distance <- amud2db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
amud2depth <- amud2db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
amud2area <- amud2db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
amud2angle <- amud2db %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
amud2distance
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 33 425. 229. 318. 51.5 313.
## 2 2 24 434. 243. 343. 57.2 345.
## 3 3 18 504. 323 435. 56.6 440.
amud2depth
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 33 38.3 10.6 22.7 7.21 21.0
## 2 2 24 56.5 12.5 30.1 12.1 28.2
## 3 3 18 98.7 27.0 59.5 23.2 51.4
amud2area
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 1 33 6317 1388 3545. 1485. 3061
## 2 2 24 11676 1998 5595. 2801. 5540.
## 3 3 18 26552 4951 12608. 6283. 11058.
amud2angle
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 33 167. 149. 158. 4.05 157.
## 2 2 24 166. 148. 157. 4.89 158.
## 3 3 18 162. 149 157. 3.78 158.
# save the results
write_csv(amud2depth, "../summary_stats/amud2depth.csv")
write_csv(amud2area, "../summary_stats/amud2area.csv")
write_csv(amud2angle, "../summary_stats/amud2angle.csv")
write_csv(amud2distance, "../summary_stats/amud2distance.csv")
# Quneitra
# Horizontal distance
qnrtdistance <- qnrtdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
qnrtdepth <- qnrtdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
qnrtarea <- qnrtdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
qnrtangle <- qnrtdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
qnrtdistance
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 21 1581 423. 947. 305. 946.
## 2 2 9 1571 595. 1071. 307. 1050
## 3 3 18 1690 710. 1072. 233. 1054.
qnrtdepth
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 21 193. 26.7 84.5 42.3 86.2
## 2 2 9 226. 35.4 106. 62.6 89.2
## 3 3 18 380. 63.4 137. 77.6 108.
qnrtarea
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 1 21 93414 8746 46069. 29058. 47146
## 2 2 9 173136 23897 71808. 54483. 45197
## 3 3 18 412713 25024 103652 93023. 69162.
qnrtangle
## # A tibble: 3 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 21 161 129. 148. 7.37 150.
## 2 2 9 166. 135. 153. 9.23 151.
## 3 3 18 160. 121. 145. 11.0 148.
# save the results
write_csv(qnrtdepth, "../summary_stats/qnrtdepth.csv")
write_csv(qnrtarea, "../summary_stats/qnrtarea.csv")
write_csv(qnrtangle, "../summary_stats/qnrtangle.csv")
write_csv(qnrtdistance, "../summary_stats/qnrtdistance.csv")
# Qafzeh
# Horizontal distance
qzdistance <- qzdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
qzdepth <- qzdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
qzarea <- qzdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
qzangle <- qzdb %>% group_by(aoi) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
qzdistance
## # A tibble: 2 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 18 1607 774. 1045. 244. 971.
## 2 2 9 1174 567 888. 208. 966.
qzdepth
## # A tibble: 2 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 18 199. 20.7 79.3 47.3 71.5
## 2 2 9 102. 9.05 47.2 34.9 39.9
qzarea
## # A tibble: 2 × 7
## aoi count max min mean sd median
## <int> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 1 18 168305 5477 46230. 40894. 32533
## 2 2 9 58111 1378 21250. 20667. 14783
qzangle
## # A tibble: 2 × 7
## aoi count max min mean sd median
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 18 155. 135. 148. 5.42 150.
## 2 2 9 157. 145. 150. 3.38 149.
# save the results
write_csv(qzdepth, "../summary_stats/qzdepth.csv")
write_csv(qzarea, "../summary_stats/qzarea.csv")
write_csv(qzangle, "../summary_stats/qzangle.csv")
write_csv(qzdistance, "../summary_stats/qzdistance.csv")
# summarise data by artefact ("site") for each measured variabled
# Horizontal distance
dbstatsdistance <- db %>% group_by(site) %>%
summarise(
count = n(),
max = max(horizontaldistance, na.rm = TRUE),
min = min (horizontaldistance, na.rm = TRUE),
mean = mean(horizontaldistance, na.rm = TRUE),
sd = sd(horizontaldistance, na.rm = TRUE),
median = median(horizontaldistance, na.rm = TRUE),
)
# Depth
dbstatsdepth <- db %>% group_by(site) %>%
summarise(
count = n(),
max = max(maximumdepth, na.rm = TRUE),
min = min (maximumdepth, na.rm = TRUE),
mean = mean(maximumdepth, na.rm = TRUE),
sd = sd(maximumdepth, na.rm = TRUE),
median = median(maximumdepth, na.rm = TRUE),
)
# Area
dbstatsarea <- db %>% group_by(site) %>%
summarise(
count = n(),
max = max(holearea, na.rm = TRUE),
min = min (holearea, na.rm = TRUE),
mean = mean(holearea, na.rm = TRUE),
sd = sd(holearea, na.rm = TRUE),
median = median(holearea, na.rm = TRUE),
)
# Angle
dbstatsangle <- db %>% group_by(site) %>%
summarise(
count = n(),
max = max(angle, na.rm = TRUE),
min = min (angle, na.rm = TRUE),
mean = mean(angle, na.rm = TRUE),
sd = sd(angle, na.rm = TRUE),
median = median(angle, na.rm = TRUE),
)
# see data summary
dbstatsdistance
## # A tibble: 5 × 7
## site count max min mean sd median
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Amud1 111 718 208. 385. 101. 370
## 2 Amud2 75 504. 229. 354. 71.3 349.
## 3 Manot 63 1163. 587. 820. 108. 823.
## 4 Qafzeh 27 1607 567 993. 241. 970.
## 5 Quneitra 48 1690 423. 1017. 281. 1008.
dbstatsdepth
## # A tibble: 5 × 7
## site count max min mean sd median
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Amud1 111 208. 5.87 38.0 24.5 31.9
## 2 Amud2 75 98.7 10.6 33.9 20.3 28.0
## 3 Manot 63 120. 25.5 62.1 18.2 57.0
## 4 Qafzeh 27 199. 9.05 68.6 45.6 59.1
## 5 Quneitra 48 380. 26.7 108. 64.6 93.0
dbstatsarea
## # A tibble: 5 × 7
## site count max min mean sd median
## <chr> <int> <int> <int> <dbl> <dbl> <dbl>
## 1 Amud1 111 30636 549 8166. 6210. 5774
## 2 Amud2 75 26552 1388 6376. 5067. 5114
## 3 Manot 63 60785 9471 27070. 10876. 25896
## 4 Qafzeh 27 168305 1378 37903. 36999. 26389
## 5 Quneitra 48 412713 8746 72489. 68399. 56725
dbstatsangle
## # A tibble: 5 × 7
## site count max min mean sd median
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Amud1 111 170. 123. 153. 8.88 155.
## 2 Amud2 75 167. 148. 157. 4.23 158.
## 3 Manot 63 160. 136 149. 5.32 150.
## 4 Qafzeh 27 157. 135. 149. 4.81 149.
## 5 Quneitra 48 166. 121. 148. 9.47 150.
# save the results
write_csv(dbstatsdepth, "../summary_stats/datastatsdepth.csv")
write_csv(dbstatsarea, "../summary_stats/datastatsarea.csv")
write_csv(dbstatsangle, "../summary_stats/datastatsangle.csv")
write_csv(dbstatsdistance, "../summary_stats/datastatsdistance.csv")
# General inventory (counts) table organised by site and aoi
inventory <- crosstable(db, c(site), by=c(aoi), total = "both")
print(inventory)
## # A tibble: 6 × 7
## .id label variable `1` `2` `3` Total
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 site site Amud1 63 (56.76%) 30 (27.03%) 18 (16.22%) 111 (34.26%)
## 2 site site Amud2 33 (44.00%) 24 (32.00%) 18 (24.00%) 75 (23.15%)
## 3 site site Manot 9 (14.29%) 12 (19.05%) 42 (66.67%) 63 (19.44%)
## 4 site site Qafzeh 18 (66.67%) 9 (33.33%) 0 (0%) 27 (8.33%)
## 5 site site Quneitra 21 (43.75%) 9 (18.75%) 18 (37.50%) 48 (14.81%)
## 6 site site Total 144 (44.44%) 84 (25.93%) 96 (29.63%) 324 (100.00%)
write_csv(inventory, "../summary_stats/inventory.csv")
Plot data
# converting to categorical
db$aoi = as.factor(db$aoi)
db$incision = as.factor(db$incision)
db$profile = as.factor(db$profile)
# Boxplot for visualise distance, depth, area and angle data organised by AOI and Site
# Angle
angle <- ggplot(data = db, aes(x = site, y = angle, fill = aoi, col = site)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Angle (º)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))
print(angle)

ggsave("../plots/angle.png")
## Saving 7 x 5 in image
# Depth
depth <- ggplot(data = db, aes(x = site, y = maximumdepth, fill = aoi, col = site)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Maximum depth (μm)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))
print(depth)

ggsave("../plots/depth.png")
## Saving 7 x 5 in image
# Distance
distance <- ggplot(data = db, aes(x = site, y = horizontaldistance, fill = aoi, col = site)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Horizontal distance (μm)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))
print(distance)

ggsave("../plots/distance.png")
## Saving 7 x 5 in image
# Area
area <- ggplot(data = db, aes(x = site, y = holearea, fill = aoi, col = site)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Area of the cut (μm2)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))
print(area)

ggsave("../plots/area.png")
## Saving 7 x 5 in image
a <- ggarrange(angle, depth, distance, area, ncol = 2, nrow = 2, common.legend = T)
print(a)

ggsave("../plots/combinedAOI.png")
## Saving 7 x 5 in image
# Boxplot for visualise distance, depth, area and angle data, together and organised by Site
# Angle
anglesite <- ggplot(data = db, aes(x = site, y = angle)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Angle (º)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1))
print(anglesite)

ggsave("../plots/anglesite.png")
## Saving 7 x 5 in image
# Depth
depthsite <- ggplot(data = db, aes(x = site, y = maximumdepth,)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Maximum depth (μm)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1))
print(depthsite)

ggsave("../plots/depthsite.png")
## Saving 7 x 5 in image
# Distance
distancesite <- ggplot(data = db, aes(x = site, y = horizontaldistance)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Horizontal distance (μm)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.8, position=position_jitter(w=0.1,h=0.1))
print(distancesite)

ggsave("../plots/distancesite.png")
## Saving 7 x 5 in image
# Area
areasite <- ggplot(data = db, aes(x = site, y = holearea)) +
geom_boxplot() +
theme_bw()+scale_fill_grey(start = 0.8, end = 1) +
labs(x = "", y = "Area of the cut (μm2)") +
scale_x_discrete(labels = c("Amud 1", "Amud 2", "Manot", "Quneitra", "Qafzeh")) +
geom_jitter(aes(colour = site), alpha=0.2, position=position_jitter(w=0.1,h=0.1))
print(areasite)

ggsave("../plots/areasite.png")
## Saving 7 x 5 in image
a <- ggarrange(anglesite, depthsite, distancesite, areasite, ncol = 2, nrow = 2, common.legend = T)
print(a)

ggsave("../plots/combinedsite.png")
## Saving 7 x 5 in image
Anova analysis
# one-way anova test the effect of raw material on damage
# The null hypothesis (H0) of the ANOVA is no difference (between artefacts) in means (for each geometric feature of the incisions), and the alternative hypothesis (Ha) is that the means are different from one another (between sites).
# check if there are differences among group means
## horizontal distance (width)
onewaymodelwidth <- aov(horizontaldistance ~ site, data = db)
summary(onewaymodelwidth)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 4 25216056 6304014 270.5 <2e-16 ***
## Residuals 319 7433069 23301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Depth
onewaymodeldepth <- aov(maximumdepth ~ site, data = db)
summary(onewaymodeldepth)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 4 208895 52224 45.39 <2e-16 ***
## Residuals 319 366993 1150
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Area of the hole
onewaymodelarea <- aov(holearea ~ site, data = db)
summary(onewaymodelarea)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 4 1.697e+11 4.242e+10 50.32 <2e-16 ***
## Residuals 319 2.690e+11 8.431e+08
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Angle
onewaymodelangle <- aov(angle ~ site, data = db)
summary(onewaymodelangle)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 4 4032 1008.1 19.41 2.6e-14 ***
## Residuals 319 16564 51.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# The p-value of the raw material variable is low (p < 0.001), so it appears that the site (at least one) has a real impact on the features.
# check normality
hist(onewaymodelwidth$residuals)

hist(onewaymodeldepth$residuals)

hist(onewaymodelarea$residuals)

hist(onewaymodelangle$residuals)

# library(car)
qqPlot(onewaymodelwidth$residuals)

## [1] 224 215
qqPlot(onewaymodeldepth$residuals)

## [1] 225 59
qqPlot(onewaymodelarea$residuals)

## [1] 225 224
qqPlot(onewaymodelangle$residuals)

## [1] 167 225
# check if model fits the assumption, the homogeneity of variances
par(mfrow=c(2,2))
plot(onewaymodelwidth)

par(mfrow=c(1,1))
par(mfrow=c(2,2))
plot(onewaymodeldepth)

par(mfrow=c(1,1))
par(mfrow=c(2,2))
plot(onewaymodelarea)

par(mfrow=c(1,1))
par(mfrow=c(2,2))
plot(onewaymodelangle)

par(mfrow=c(1,1))
# The model fits the assumption of homoscedasticity
# ANOVA tells us if there are differences among group means (sites and features), but not where those differences are. To find out which groups are statistically different from one another, you perform a Tukey’s Honestly Significant Difference Tukey’s HSD post-hoc test for pairwise comparisons.
# check which groups are statistically different from one another
TukeyHSD(onewaymodelwidth, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = horizontaldistance ~ site, data = db)
##
## $site
## diff lwr upr p adj
## Amud2-Amud1 -31.13481 -93.73121 31.46158 0.6508190
## Manot-Amud1 435.00665 368.94828 501.06502 0.0000000
## Qafzeh-Amud1 607.55586 517.69281 697.41890 0.0000000
## Quneitra-Amud1 631.97252 559.62870 704.31634 0.0000000
## Manot-Amud2 466.14146 394.57268 537.71024 0.0000000
## Qafzeh-Amud2 638.69067 544.70264 732.67870 0.0000000
## Quneitra-Amud2 663.10733 585.69926 740.51541 0.0000000
## Qafzeh-Manot 172.54921 76.22087 268.87755 0.0000141
## Quneitra-Manot 196.96587 116.73240 277.19934 0.0000000
## Quneitra-Qafzeh 24.41667 -76.32592 125.15926 0.9637087
plot(TukeyHSD(onewaymodelwidth, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodeldepth, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = maximumdepth ~ site, data = db)
##
## $site
## diff lwr upr p adj
## Amud2-Amud1 -4.119425 -18.028363 9.789513 0.9266782
## Manot-Amud1 24.027616 9.349427 38.705805 0.0000966
## Qafzeh-Amud1 30.536664 10.569067 50.504260 0.0003391
## Quneitra-Amud1 70.100691 54.025873 86.175510 0.0000000
## Manot-Amud2 28.147041 12.244437 44.049646 0.0000185
## Qafzeh-Amud2 34.656089 13.771920 55.540257 0.0000736
## Quneitra-Amud2 74.220117 57.020019 91.420215 0.0000000
## Qafzeh-Manot 6.509048 -14.895138 27.913233 0.9198008
## Quneitra-Manot 46.073075 28.245175 63.900976 0.0000000
## Quneitra-Qafzeh 39.564028 17.178994 61.949061 0.0000191
plot(TukeyHSD(onewaymodeldepth, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodelarea, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = holearea ~ site, data = db)
##
## $site
## diff lwr upr p adj
## Amud2-Amud1 -1790.275 -13697.294 10116.74 0.9939065
## Manot-Amud1 18903.955 6338.403 31469.51 0.0004484
## Qafzeh-Amud1 29737.055 12643.403 46830.71 0.0000272
## Quneitra-Amud1 64322.435 50561.272 78083.60 0.0000000
## Manot-Amud2 20694.230 7080.494 34307.97 0.0003762
## Qafzeh-Amud2 31527.330 13649.029 49405.63 0.0000201
## Quneitra-Amud2 66112.710 51388.229 80837.19 0.0000000
## Qafzeh-Manot 10833.101 -7490.372 29156.57 0.4844187
## Quneitra-Manot 45418.480 30156.556 60680.40 0.0000000
## Quneitra-Qafzeh 34585.380 15422.233 53748.53 0.0000118
plot(TukeyHSD(onewaymodelarea, conf.level=.95), las = 1, cex.axis=0.5)

TukeyHSD(onewaymodelangle, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = angle ~ site, data = db)
##
## $site
## diff lwr upr p adj
## Amud2-Amud1 3.9949550 1.039992 6.9499179 0.0022690
## Manot-Amud1 -4.4926641 -7.611055 -1.3742732 0.0009003
## Qafzeh-Amud1 -4.7894895 -9.031618 -0.5473606 0.0179890
## Quneitra-Amud1 -5.3304617 -8.745567 -1.9153561 0.0002366
## Manot-Amud2 -8.4876190 -11.866138 -5.1091004 0.0000000
## Qafzeh-Amud2 -8.7844444 -13.221300 -4.3475892 0.0000011
## Quneitra-Amud2 -9.3254167 -12.979589 -5.6712446 0.0000000
## Qafzeh-Manot -0.2968254 -4.844159 4.2505079 0.9997691
## Quneitra-Manot -0.8377976 -4.625347 2.9497516 0.9739556
## Quneitra-Qafzeh -0.5409722 -5.296687 4.2147428 0.9979346
plot(TukeyHSD(onewaymodelangle, conf.level=.95), las = 1, cex.axis=0.5)

# Some observations from the post-hoc test results:
# All raw materials are significantly different from each other. All p-value (here represented as p adj) are smaller than 0.05.
End of Script